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ABSTRACT 


A new technique for estimation of the spectrum of an 
unknown Signal is investigated. The technique combines a 
radiometric method of signal detection (modulating the re- 
ceiver noise with the assumed signal by alternately connect- 
ing and disconnecting the input transducer) with a Fast 
Fourier Transform processor for spectral analysis. After 
obtaining averaged spectral estimates with the transducer 
(antenna) connected and disconnected the difference is formed 
and statistical analyses applied to determine statistically 
Significant terms. The significant components comprise the 
estimate of the spectrum of the signal. 

Analytical results are derived for sinusoidal signals and 
uncorrelated, Gaussian noise. Experimental results for both 
narrow and wideband signals and additive Gaussian noise are 
presented. Reasonable spectral estimates of a wideband signal 
are obtained with 16 averages when the composite signal-to- 


noise ratio is O dB. 
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I. INTRODUCTION 


Wideband noise-like signals occur in various types of 
radio communications and radio ranging systems. In this re- 
port, we present a method of signal processing which provides 
a measure of the spectrum of a wideband signal. This method 
combines the technique of radiometry which is used to detect 
weak radio signals with the Fourier transform which defines 
the spectrum of an electrical signal. 

For more than twenty years, radiometers have been used in 
radio astronomy and elsewhere to detect low power signals 
(weak stars). In a radiometer the receiving element (antenna) 
1s alternately connected to and disconnected from the receiver 
as shown in Fig. 1. The received voltage is therefore alter- 
nately weak signal s(t) plus noise n(t) and noise only. Hence, 
the receiver noise serves as a "carrier" which is then ampli- 
tude modulated by the weak signal by virtue of the switching 
action. Typically, a square-law device is used to demodulate 
the carrier. The square-law detector makes the system a power 
detector. Integration provides energy detection capability. 
The output from the energy detector alternately represents 
the energy level with no signal present (antenna disconnected) 
and the energy level when it is possible that a signal is 
present (antenna connected). The two values are compared 
uSing appropriate statistical analysis and a decision is 


made concerning the presence or absence of a signal. [Ref. 1] 
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So, the use of the receiver's own noise as a reference 
makes possible the detection of weak signals when the system 
noise 1s predominantly receiver noise (thermal noise-limited 
receivers). [In general, systems are thermal-noise limited 
if the frequencies of interest exceed 20 MHz.] Our interest 
fom che Measurement of the spectra of weak signals. Re- 
placing the square-law detector of Fig. 1 with a Fourier 
transform processor results in the system of Fig. 2. (Ina 
practical system, the receiving subsystem of Fig. 2 includes 
a mixer so that the FFT processor operates on a lowpass or 
baseband portion of the spectrum.) 

In this report, the processor of Fig. 2 includes a sampler 
which then allows use of the Fast Fourier Transform (FFT) 
algorithm to generate a discrete form of the spectrum of the 
processor input voltage. The FFT processor output iS manipu- 
lated to produce an array of numbers proportional to the power 
at discrete frequencies, (the discrete power spectrum of the 
input). An average over M blocks of input data is formed which 
provides an output proportional to the energy of the input. 

Because of the switching action at the system input, the 
FFT processor output alternately represents the spectrum 
Gy (£) of the noise and the spectrum Go (£) of the combination 
of signal and noise. For the case of additive noise uncorre- 
lated with the signal, we can write the signal spectrum Go (f) 
as Go (f) = Gy (ft). Therefore, by subtracting the averaged 


spectra obtained with the FFT processor, a measure of the 
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spectrum of the signal is formed when a signal is present. 

A near-zero difference indicates the absence of a signal. 
Fig. 3a shows the spectrum Gy (f) @f the me@1se only, Fig. 3b 
shows the spectrum G. (f) of the combined signal and noise for 
a square wave signal, and Fig. 3c shows the difference Go (f). 
Statistical analyses relate the number of blocks of data used 
in averaging to false alarm rate and detection probabilities. 
In this manner, the significant values of the difference 
(signal) spectrum are determined. Fig. 3c shows that the 
seventh and higher order harmonics of the square wave are 
missing from the spectral estimate. This is because the 
energy content of the higher order harmonics cannot be dis- 
cerned from the energy content of the noise. Fig. 3c also 
shows two components which are not part of the spectrum of 
the signal (one at approximately 30 and another at approxi- 
mately 150). These are caused by false alarms in the statis- 
tical analysis. 

In this report, we identify this method with the acronym 
FABLUS (Fourier Analysis of Broadband, Low-level Signals 
Using Switching). 

The FABLUS method of estimating signal spectra is useful 
when the receiver noise limits system performance. Noise which 
appears at the antenna (transducer) terminals is treated as 
Signal. For accurate estimates of signal spectrum this method 
also requires the noise to be stationary over the interval of 


time required to obtain all the data used in calculating the 
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averages of the spectra of noise and combination of signal and 
noise. When these two conditions are satisfied, there is no 
theoretical limit on the level of a signal which can be de- 
tected and the accuracy with which the spectrum can be mea- 
sured. By averaging over a large number of blocks of data 
(many ensembles or many cycles of the switch), the spectrum 

of a comparatively very weak signal can be accurately measured 
by this method. 

Table 1 is an example of the results obtained using this 
method. Shown is the probability of detecting a signal com- 
ponent as a function of the ensemble size M and the ratio of 
the average signal amplitude to the average of the noise 
amplitude (frequency interval signal-to-noise ratio) when the 
false alarm rate is 2.5 percent. 

The concept of signal-to-noise ratio requires definition 
when using the discrete Fourier transform. For example, if 
one watt of white, Gaussian noise (bandlimited to the pro- 
cessor bandwidth) is applied to an FFT processor which uses 
a 1024 sample block length the output will contain 1/512 watts 
in each frequency interval. If a one watt sinusoidal signal, 
centered on one of the frequency intervals, is applied to 
the processor the output at the interval under consideration 
will contain one watt of power. The signal-to-noise ratio at 
this interval is 1/(1/512) or 27.1 dB. However, the signal- 
to-noise ratio at the input to the processor is 1/l or 0O dB. 


So, in evaluating the performance of a detection method, care 
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must be used in defining the signal-to-noise ratio. All of 
the signal-to-noise ratios contained in this report are fre- 
quency interval (bin) signal-to-noise ratios based on a 
block length of 1024. 

Chapter II of this report contains a detailed discussion 
of FABLUS including analysis of the errors associated with 
application of the FFT and the operations necessary to pro- 
duce the arrays of numbers representing the energy content 
of noise only and combination of signal and noise. Chapter 
III discusses the statistical techniques investigated to 
determine the significance of the spectral components. Experi- 
mental results obtained using FABLUS are contained in Chapter 
IV. Chapter V presents conclusions. Derivations of analyti- 
cal results are contained in Appendices A, B and c. A list 


of references is provided. 
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ti. Sob PrABLUS PROCESS 


Pee © ROCEDURE 

Since the FFT processor is a digital device, the analog 
input voltage must be sampled and quantized. Fig. 4 shows an 
analog waveform x(t) and the sampling operation. (The errors 
caused by quantization are not addressed in this report.) 
The sampling interval T and the number N of sample points 
in a data record are values which must be chosen carefully 
if the FFT processor output is to represent the spectrum of 
the input. Appendix A contains a brief review of Fourier 
transform theory and the Discrete Fourier Transform (DFT). 

The output of the FFT processor is shown in Fig. 5 where 
X(m£ >) represents the discrete Fourier transform of x(t); 
mis the set of integers. The output consists of N complex 
numbers; however, only the first N/2 are unique. The spectra 
is "folded" about the N/2 value. Therefore, the usable output 
consists of N/2 complex numbers representing the inphase and 
quadrature components of the transform. These components 
are the transform values at discrete frequencies at intervals 
of 1/NT Hz in the range 0 to 1/2T Hz. We refer to an interval 
of 1/NT Hz as a frequency bin. We denote the output numbers 


by G, when the combination of signal and noise is the input 


C 
and by Gy when noise only is the input. These numbers vary 
from one block of data to the next because of the random 


nature of the input (random phase) relative to the switching 
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and data truncation which is part of the FFT processor opera- 
tion. To remove this randomness the squared magnitude of the 
transform (Re ( \4 + Im( )4) 1s computed. The squared magni- 
tudes are insensitive to phase variation and are proportional 
to the power contained in the input time record at each fre- 
quency bin. Fig. 6 shows the FFT output magnitude squared. 
Therefore, the first processing step is to convert the FFT 
output to magnitudes squared (16,14 and 1G 17) at each bin. 
Averaging of M output arrays smooths the magnitude squared 
values and provides values proportional to energy for each 
bin. To perform the averaging, the sum of the magnitudes 
Squared must be computed and stored after each block of input 
data has been processed. In addition to storing the sum of 
the magnitudes squared, the sum of the squares of the magni- 
“lle is also stored to allow 


of the signal plus noise data 


tudes squared (1|G,]*1° and [|G 


computation of the variance one 


and of the variance =e of the noise only data. Since the 


Switching provides alternately noise only and a combination 
of signal and noise, a total of 2M blocks of data must be 
processed before computation of means (\G,|* and 1G |) and 


sample variances (s : and s.°) can be done. In the remainder 


& 
of this report M will be referred to as ensemble size and 
will refer to the taking of a total of 2M samples of input 
voltage. These four arrays are the numerical output from 
which the estimate of the spectrum of the signal is formed. 
ie 2 7 


1 So and Sy) 








Appendix B contains derivations of icale.. |G, 


for the case of zero mean, white, Gaussian noise. 
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B. PROCESSING ERRORS 

The output of the FFT processor contains the information 
necessary to estimate Go(f), the spectrum of s(t). Variations 
in the arrays ex and Go due to random phase have been elim- 
inated by using the transform magnitude squared. However, 
there are still four potential sources of error which could 
cause this estimate of G. (f) to differ from its true value. 
Three of these error sources are sampling of the input wave- 
form, truncation of the input data record and the discrete 
nature of the output spectrum. These three are all associated 
with the FFT process. The fourth potential source of error 
1s associated with the switching action necessary to produce 
the alternating inputs of noise only and combination of signal 
and noise. 

Appendix A contains a discussion of the errors associated 
with estimating a continuous spectrum using an FFT processor. 
As derived there, the output of the FFT processor, with input 


x(t) can be viewed as composed of three terms: 
= 1 —_ “= 
X (mE) ) = xX (mf) Al (m£)) Tr (mf) 


In this equation X (mf) is the estimate of X(f), Al (mf) 
is the error caused by aliasing in the sampling of x(t), 
Tx (mf) is the error caused by truncation of the input data 
and X' (mf) is the sampled replica of X(f), the continuous 
spectrum; fo is the frequency interval = 1/NT. If the input 


to the FFT processor is noise n(t) only, the output will be 


29 








Gy (mE 9) = G' (mEq) - Al, (mf) - Try (mf) 


When the input is signal plus noise, s(t) + n(t), the pro- 


cessor output becomes: 
— t a =_ 
G, (mf )) = G, (mf) Al, (mf,) Tr. (mf) 


The errors due to aliasing in the signal plus noise output 
can be considered as composed of components associated with 
the noise (Al) and components associated with the signal 
(Al). If the signal and noise are uncorrelated, these two 
components will be independent and Al. (mf) can be written 
as Al. (mf£)) + Al, (mf) - Errors associated with the trunca- 
tion of the input for signal plus noise conditions can also 
be written as two components Tr, (mf 9) = Tr. (mf)) + Try (m£ 9) 


because of the linearity of the Fourier transform. Form 


R(mf£) Go (mt) = Gy (mt) 


- t is aa t 
Ge. (mf) Al. (mf) Tr. (mf£,) - Gn (mE) ) 


= Al, (mf) = Try (mf 9) = G.' (mf) Al, (mf ,) 
+ Tr), (m£ 9) 
= Go (mf) - Al. (mf)) = Tr. (mf) 


This result shows that the errors caused by aliasing and 
truncation of the noise only waveform are eliminated when the 
difference is formed. Therefore, the computed difference 


represents the spectrum of the signal as if it had been applied 
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to the FFT processor without the noise. Since R(mf£) SCen— 
tains these erroneous frequency components which are not 
eliminated by the post transform processing, care must be 
exercized to minimize their effects. Careful choice of T, 
the sampling period, and N, the number of samples per block, 
Will accomplish this. 

One of the major features of this detection scheme is 
the switching of the antenna into and out of the receiver 
system input. This switching action can be modeled as a 
Squarewave alternating between 0 and +1. Call this waveform 
j(t). With s(t) the signal and n(t) the noise the FFT pro- 


cessor input x(t) is: 
MG) Ste} git) + nt) 


The multiplication of s(t) by j3(t) produces discontinuities 
in x(t). Multiplication in the time domain corresponds to 
convolution in the frequency domain. Denoting the convolution 
operation by *, we then have that the spectrum X(f) of x(t) 


is: 
ea) a eee on) *) (tf) + N(f£) 


The convolution of J(f£) with S(f) produces frequency components 
in X(f) which are not in S(f). However, these spurious com- 
ponents can be eliminated by blanking the FFT processor input 
during the transitions of j(t). This will require that the 


record length of the processor input be slightly less than 
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Piewiart pertea of F(t) and that the FFT input blanking pulse 
and j(t) be synchronized to allow the transients to decay 
sufficiently before the data is applied to the FFT processor. 
That is, the effects of switching are eliminated by selecting 
a complete block of data to be processed during the time 
interval the switch is closed. 

In summary, the output of the FFT processor 1S an array 
of numbers representing an estimate of the spectrum of the 
processor input. This estimate is influenced by the phase 
variations of the data blocks relative to the signal, and 
by sampling, truncation and switching. The random phase 
variations are removed by computing the magnitude squared. 
The errors associated with the switching are removed by 
synchronization of the data taking and switching operations. 
Post transform processing (differencing) eliminates errors 
caused by sampling and truncation of the noise. The errors 
remaining in the post transform processed data are those 
associated with sampling and truncation of the signal s(t). 
These are errors in the FFT process and are not caused by 
the switching action. 

Fig. 7 1S a more detailed block diagram of the FABLUS 


system. 
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IIIf. STATISTICAL METHODS 


The output of the FFT processor is four arrays of numbers 
representing the sum of the squared magnitudes and the sum 
of the squared magnitudes squared (fourth power) at each fre- 
quency bin for noise only and signal plus noise inputs. The 
next step in the detection process is to manipulate these 
arrays to perform a statistical evaluation at each bin. The 
evaluation determines which bins contain statistically signi- 
ficant components. Those bins with significant components 
and the magnitudes of those components comprise the spectral 
estimate of s(t). 

There are at least two general courses of action available 
to perform the statistical analysis. First, it can be assumed 
that a signal is present and then evaluate the arrays to 
confirm or deny this hypothesis. Alternatively, it can be 
assumed that no signal is present and then evaluate the arrays 
to confirm or deny this hypothesis. These two approaches are 
not equivalent. 

Hypothesizing that a signal is present requires some 
assumption about the level of that signal. This assumed level 
is the smallest which can be detected with the desired degree 
of confidence. Then, the statistical evaluation determines 
with l-a confidence that a signal of the assumed level or 
higher exists. This procedure sets a detection rate and a 


rate of missed detections. 
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Hypothesizing that no signal is present requires no 
assumptions about the signal level. The statistical evalua- 
tion determines with l-a confidence that no signal is present. 
This method sets the rate at which correct no signal decisions 
are made and the false alarm rate. The relation between these 
hypotheses, and the detection and false alarm rates are de- 
Beeecca in Fig. 8. 

The no signal hypothesis requires no assumptions about 
Signal level and setting of the false alarm rate. Unfortu- 
nately, the signal detection rate cannot be determined with- 
out signal level assumptions. 

To minimize the number of assumptions required, the hypo- 
thesis no signal present is used in the signal detection 
methods employed here. This is in agreement with general 
practice where control of false alarm rate is usually pre- 
ferred to control of signal detection rate. Various methods 
are investigated. The two most useful methods are first des- 


Cribed in detail in the following. 


A. METHOD I--DIFFERENCE OF TWO MEANS 
This method requires fewest assumptions. The only assump- 
tion is that M (the number of samples used in acquisition of 
the arrays stored by the FFT processor) be large enough so 
ig 


the sample means (|G. a) are approximately Gaussian 


distributed and the sample variances (eee and s.°) are approxl- 


and |G 


mately the population variances (cs and 00°) - Sample sizes 
as small as 10 or even less can meet the first requirement 


to reasonable accuracy. It is more difficult to predict the 
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number of samples required to meet the second condition. 
Values of 30 or 40 are cited in some texts [Ref. 2,3]. How- 
ever, Since the population variances are unknown, no other 
reasonable action iS available and this substitution is used 
for all values of M. 


Under the assumptions listed above, the random variables 
\g 








|G 


actual mean of the noise population) and v (the actual mean 


and isl? can be considered Gaussian with means u (the 


of the signal plus noise population) and variances on /M and 





on /M respectively. Form the random variable R= IGe|* = 





|G which is Gaussian distributed with mean (v-u) and 


2 


n | 


variance (Ca, + O°) /M) « As indicated above the population 


variances must be replaced with the sample variances to yield: 


2 


Random Variable R: Normal ((v-u) , (S, +s.) /M) 


Without loss of generality we can assume that no signal is 
present (v = u) and make the null hypothesis HO: u =v. 
The method proceeds as follows: 


For each frequency bin 1: 





“2 a2 2 Z 
a.) Compute Gu ae Go i, Sy and Sai 
b.) Compute R. = |G |i? - |G ie and 
| : iH GC N 
a 
S1 Ni 


So; = (Sa;° + Sys )/M 


Si 


| c.) Compute Cv = zs 
where zy = the standardized normal variable such 


that the area under this density function from 


ZO to ~ equals a. 
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| 


Ge Lt Rs < Cv accept HO and make decision DO = no 


Signidailwwat Din 1. 


ea) Lt R. >= Cv reject HO and make decision Dl = signal 


present at bin 1. 


f.) If HO is accepted set R; = 0; spectral estimate 


does not have a component at ith bin. 


g-) If HO is rejected, then R, 1s the estimate of the 
magnitude of the spectral component of the signal 
aimee 2. 
Fig. 9 is a block diagram of this method. 
If no signal is present then the random variable R, should 
have zero mean with variance 254° /M Since the variance is 
the same for both sets of M samples. The probability density 
fmec ion for R, conditioned on no signal is: 
pee o— vie = /Mi/ (2s, Vm) exp (-Mr*/4s,*) 
where r is a value of the random variable R.. The critical 
value is Zz V2s,/ 7M. The region to the left of the critical 
value represents those values of R for which HO is accepted; 
to the right, HO is rejected. From this it is obvious that 
tame choice of zy sets the false alarm rate. 
If a Signal is assumed to be present, then v-u = d and 
2 


Ry has mean d and variance (Sa; +s. 67) /M. The probability 


density function for Ry conditioned on a signal present is: 


— / 2 2 2 ; 2 2 
poeekv-u=d) = /M/ (27(S,,° +8), ) )exp(-M(r-d) “/2(s,,°+S,5°)) 


a 
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The critical value is not the same as in the no signal case 
because the variance of R is greater due to the signal. The 
probability of detection will be a function of d and Sc: 
The value of ZO is set by the false alarm rate. Fig. 10 is 

a graph showing both the no signal and signal present density 
functions. 

To evaluate the analytical performance of this spectral 
estimation method, we apply the results from Appendix B to 
the procedure described above. The probability of detection 
for Method I is P(r >Cv) where Cv is the critical value. 


This is evaluated in Appendix C for the case of sinusoidal 


Signal and zero mean, white, Gaussian noise to give: 
P(r >Cv) = Q(z, - (A*N*/M) /4V(N*ACNN*/2 + 2N* (n*) *))) 


In this equation A is the peak amplitude of the sinusoidal 
Signal, n2 is the mean square value of the noise process, M 

is the ensemble size, N is the number of points in the input 
data block and zy 1s, as before, the value of the argument of 
the normal probability density function such that the area 
under the function from z, to » is a. The function Q(x) is 
defined as the area under the standard normal density function 
from x to t~. Fig. 11 is a graph of probability of detection 
versus bin signal-to-noise ratio. From Fig. 1l we see that at 
0 dB bin signal-to-noise ratio, the probability of detection 


increases by 22 percent (from 30 to 52) when the number M of 


records averaged is doubled (from 8 to 16). 
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B. METHOD II--MEANS ONLY 

Method I requires that the mean and variance of both the 
Signal plus noise and the noise only ensembles be computed. 
These computations can be time consuming. Method II was 
developed to reduce the data and computation requirements. 
For this method only the mean values of the signal plus 
noise and noise only ensembles are required. However, the 
assumptions required are more restrictive than Method I. 
In addition to the assumption that the sample means are Gaussian 
distributed, it is assumed that a linear relationship exists 


between the standard deviation (Ss) of the noise only and the 








mean afer |) of the noise only, and that the variance of both 
ensembles are equal. That is, it is assumed sy = Bees 
and o 2 Oo 2 where B is a constant. 


N S 


It is difficult to identify cases where this linear rela- 
tionship exists unless analytical methods are employed. 
Appendix A shows that for Gaussian, white noise these assump- 
tions are valid with B = 1. The equal variance assumption 


is in general not valid; however, for small signal levels 











(a —< -), the error introduced is small. 

Under these assumptions, the random variables IG |? and 
je,|* are Gaussian distributed with means u and v respectively 
and common variance on. Then the random variable Rs = Ico? - 





1" is Gaussian with mean (v-u) and variance 20,,°/M. With 


the above assumptions we have: 


Random Variable R: Normal((v -u), 2B* (|G |7)7/m) 
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As in Method I we assume that u = v and the null hypothesis 
HO: u = v. 


For each frequency bin i: 


ae ae 
a.) Compute Go 4 and ere: 
b.) Compute R. = |G 1a? - |G 1a? and 
an Cc N 
2,55 
Sa, = VSB|G. |i / VM 
eee Compute Critical value = Cy = Zz Soi 


ee) «61 E R; < Cv accept HO and make decision DO 


no 


Signal at bin i. 


e.) If R; > Cv reject HO and make decision Dl = sig- 


nal present at bin i. 


f.) If HO is accepted set R, = 0; spectral estimate 


does not have a component at ith bin. 


g-) If HO is rejected, then R, is the estimate of the 





magnitude of the spectral component of the signal 
| aceoin 1. 
Fig. 12 is a block diagram of the method. 

Method II is basically the same as Method I. The only 
substantial difference is the method of establishing the 
variances and therefore the critical value. 

Assuming no signal is present the decision variable R is 
GausSian distributed with zero mean and variance 25° /M. 
Therefore the probability density function of R conditioned 


on no signal is: 
= Z 2 
P(rju=v) = VM/(2s).¥n) exp(-r M/4s.°) 


40 










3 
| in SIGHAL 
. PRESENT 


IF O36 
+ 0 
SIGNAL NOT 


PRESET 
iF be 


J 


wig tes ee lock diawram of detection Method I1. 


41 





substituting BIG, |° for s. 





Pac) /Mi/ (2B|G. | “Vm) exp (-Mx*/(4B°{ |G, |“)*) 


The critical value will be: 





- p 
Cv = Zz V2B(G, | /VvM 


Like Method I the area under the density function to the left 
of the critical value represents the probability of accepting 
HO (correct no signal decision). The area to the right repre- 
sents the probability of incorrectly rejecting HO (false alarm). 
If B is known exactly, then Za will set the false alarm rate. 
If a signal is present the decision variable R. is Gaus- 
Slan distributed with mean (v-u). However, the variance 1s 
the same as in the no signal case. Therefore, the density 


function conditioned on signal present will be: 





P(rjv-u=d) = // (2B |G, |“Vm) exp (-M(x-d) */4B°{ |G. |°)*) 


The area under the curve to the left of the critical value 
represents the probability of accepting HO (missed detection). 
The area to the right is the probability of correcting de- 


tecting the signal. In this case the probability of detection 
re 





will be a function of d,B and |G Fig. 13 shows the 


N 


density functions under signal and no signal conditions. 


C. OTHER METHODS 
In addition to the methods described above, we also in- 
vestigated two methods which require the assumption that the 
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population being sampled is Gaussian distributed. The two 
methods are Smith-Satterthwaite test for the difference of 
two means using Student t statistics and the ratio of two 
variances using F distribution statistics. The results for 
both these methods are significantly poorer than those of 
either Method I or II. 
1. Smith-Satterthwaite Test--Difference of Means 

This test suffered from two identified shortcomings. 
Unlike the two previous methods where the accuracy of the 
assumptions increased with increasing M, the validity of 
assuming the population being sampled is Gaussian distributed 
is not a function of M. Therefore, errors caused by using 
this assumption persist for any ensemble size. In addition, 
because the critical value is determined by calculating the 
degrees of freedom from the variance values, it is difficult 
to get consistent false alarm performance [Ref. 4]. 

2. F Distribution--Ratio of Variances 

The ratio of two variances test also is afflicted 
with the errors caused by assuming the sampled population 
is Gaussian distributed. In addition the variability of the 
ratio of two variances is quite large at the relatively small 
ensemble sizes of interest here. This coupled with the 
rather large right hand tail of the F distribution causes 
the performance to be substantially worse than any of those 


already discussed [Ref. 5]. 


Ad 





D. REFERENCE CALCULATION 
One of the serious sources of error in Methods I and II 
is the variability of the critical value Cv. As long as the 


process is stationary the critical value should be a constant: 


ae 
Since the population variances are unknown, we use S, as an 
estimate of Tas The changing value of Cv will cause errors 


in detection in addition to those due solely to the random- 
ness associated with the decision variable R, . 

To gain insight into the amount of error caused by the 
changes in the critical value, this reference method was 
developed. To establish a constant critical value, the values 
of the variance of the noise and signal plus noise ensembles 
must be known exactly. To obtain a good estimate of these 


variances the statistics of both the noise and signal plus 





noise ensembles were collected during the testing of Methods 
I and II. This information was then evaluated to establish 


the values of Bee and Onn based on all the data taken during 
the tests. These numbers were used to calculate ee The 


computation proceeds as follows: 
At each frequency bin i: 
a.) Compute |G, |i and Ie, |i? 


2 


2 
b.) Compute R, = IG. |i 


aCe 


c.) Read in o, and compute Zo 


€ e 
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a) Peake < Zo 
ale 


oC then accept HO; no signal at bin i 


e.) If R; > Zid reject HO; signal present at bin i 


C 
f.) If HO accepted set R, = 0; spectral estimate has 


no component at bin i. 


Because this method eliminates the errors caused by 
variability in the critical value, the detection results 
Should be the best obtainable for a given signal level. 
Thus the detection results obtained with this method are 
used as a reference to compare the performance of Methods I 
amet. Fig. 14 is a block diagram of the computation pro- 


cess used to obtain the reference results. 


E. SIGNAL-TO-NOISE RATIOS 

The bin signal-to-noise ratio is the ratio of the average 
Signal magnitude squared to the average noise magnitude 
squared Mee-7ic,1°) for a particular frequency bin. For 
the case considered above the signal-to-noise ratio is, from 


Appendix A, 


SNRL aAN*/4Nn* 


ACN/4n* 


This is a function of N, the input data block length, when 
using the FFT algorithm. 
Alternatively an expression for the composite signal-to- 


noise ratio can be developed by solving 


IG{* = Nn* and IG 
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fey a and a and forming: 


SNR2 = ye 





2 2 
2|G,|“/N|G,.| 
In this case the factor N is in the denominator whereas be- 
fore it was in the numerator. For N = 1024 the following 
relationship exists between these alternative representations 


when SNR1L and SNR2 are in decibels and where 10 log N/2 = 


10 log 512 = 27.1. 
SNR2 = SNRI1 - 27.1 GB 


The expression shown for SNR1 is used throughout this 
report since the methods employed make spectral estimates on 
a bin-by~bin basis. Consequently, the noise power (or energy) 
at a particular frequency bin influences only the spectral 
estimate at that bin. Furthermore, in the general case of 
wideband, noise-like signals and non-Gaussian noise, SNR2 is 

2 


not calculable because formulae for A’ and n* as functions 


of the average magnitudes squared are not known. 
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IV. EXPERIMENTAL RESULTS 


A. NARROWBAND RESULTS 

To demonstrate the capabilities of this method, a hard- 
ware evaluation was conducted using a signal generator, noise 
source, spectrum analyzer and minicomputer. A square wave 
Signal was used to allow simultaneous processing of multiple 
Signal levels. The noise source produced white, Gaussian 
noise. The spectrum analyzer did the sampling and obtained 
the FFT. The minicomputer provided the switching control, 
storage and computational capacity necessary to carry out the 
statistical analyses. 

A total of 1000 trials was conducted with ensemble size 
M of eight. After each set of data was collected and the 
various sums computed, the values at the frequency bins con- 
taining the primary and third, fifth, seventh, ninth and 
eleventh harmonic frequencies of the input square wave were 
stored as well as data from five bins which contained no sig- 
nal. After the 1000 trials were completed the data collected 
was analyzed using the previously described Methods I and II. 
A preliminary analysis was conducted using z= 1.96 for Method 
I and zy = 1.96 and B = 1 for Method II (2.5 percent false 
alarm rate). Because of the variation in the validity of 
the assumptions used, the desired false alarm rate was not 
achieved. The value of Zz was then adjusted separately for 


each method until a 2.5 percent false alarm rate was achieved. 
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Next the data was reprocessed to determine the average values 
and variances at each bin over the 1000 trials. MThese aver- 
ages and variances were then used to process the data to ob- 
tain the reference results. Table 2 shows the detection 
results for these trials. 

The data obtained for the 1000 trials with M = 8 was then 
reprocessed to yield data for 500 trials with ensemble size 
16. This data was then processed using Methods I, II and 
the reference procedure. These results are presented in 
Table 3. 

The results for both the M = 8 and M = 16 trials are 
plotted in Fig. 15 along with the analytical probability 
of detection curves. 

A comparison of the analytical and experimental results 
reveals that: 

1. The analytical and experimental results for 
Method I agree very closely. 

2. The reference results are substantially better 
than either Method I or II. 

3. Method I is slightly better than Method II. 

The very close agreement between the analytical and ex- 
perimental results for Method I indicates that the assumption 
used (sample mean is Gaussian distributed) is generally valid 
for the case of white, Gaussian noise and sinusoidal signals. 
The results are not expected to change for other signal for- 


mats, Since any waveform encountered in practice can be 
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TABLE 2. Probability of detection (2%. Ensemble size 3, 
false alarm rate 2.5%. 





PRG@phetEITY OF DETECTION ¢x%) 
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Dimeee S. Prebability at detection (4). Ensemble size 15, 
false alarm rate 2.5%. 
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expressed as a sum of sinusoids. The validity of the assump- 
tion for other types of noise is open to question. However, 
the ability of the method to detect the signals will not be 
greatly changed by the presence of non-Gaussian noise if the 
ensemble size is large enough to allow the Gaussian dis- 
tributed sample mean assumption. 

As expected the reference results show greater detection 
ability at all signal levels. They are included here only 
to indicate an upper limit of detection under these conditions. 
It may be difficult to implement such a method since it re- 
quires prior knowledge of the statistics of the processes 
involved. 

A comparison of the probabilistic characteristics of 
Methods I and II (see Figs. 10 and 13) indicates that Method 
II should have higher probability of detection than Method I 
because the critical value will be less under otherwise 
identical conditions. However, the experimental results show 
precisely the opposite. The smaller critical value obtained 
using Method II is due to the use of the noise only data to 
compute the variance of R. The resulting smaller critical 
value does increase the probability of detection, but it also 
increases the false alarm rate. When the value of zy 1s ad- 
justed to regain the desired false alarm rate, the probability 
of detection decreases to the levels shown in Tables 2 and 3. 
The increase in the false alarm rate is not predicted by the 


analytical analysis but is due to additional noise entering 
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the experimental system along with the signal. This type 

of error can be encountered in any method where the critical 
value is calculated from noise only data. The use of both 
Signal plus noise and noise only data to calculate the criti- 
cal value, as in Method I, provides a limited amount of 
protection against this type of error. 

The marginal increase in detection capability which 
Method I shows over Method II is probably due to the greater 
amount of information used in the Method I process. Since 
Method I uses both the means and variances to make spectral 
estimates, the extremes in the means may be counterbalanced 
to some extent by the variances. This might produce a smaller 
overall fluctuation in the decision variables. Method II 
uses the means to set the critical value and to compute the 
value of the decision variable. Therefore, this process 
would be more sensitive to the variations of the means. 

However, the slightly increased detection capability of 
Method I does not occur without cost. Implementation of 
Method I requires additional computation time since the arrays 
of squared magnitudes squared must be computed and requires 
twice the storage capacity since four arrays must be stored 
instead of two. 

A graph of ensemble size versus signal-to-noise ratio for 
50 and 95 percent probability of detection is included as 
Htayet6. Additionally, Fig. 17 is a graph of probability of 


detection versus bin signal-to-noise ratio for 2.5, 5 and 10 


So 
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percent false alarm rates. Figs. 16 and 17, which are ob- 
tained from the equations derived in Appendix B, indicate 


the trade-offs involved using these detection methods. 


B. WIDEBAND RESULTS 

The FABLUS process makes decisions on a bin-by-bin basis, 
and so it can be applied to wideband signals having continu- 
ous spectra. To test this capability, a pseudo-random se- 
quence of length 1,048,575 was used as a signal. Because a 
noise-like signal is used, the signal spectrum can only be 
specified on an average or mean basis. Fig. 18a shows the 
spectrum of the pseudo-random sequence based on an average 
over 500 data blocks. Fig. 18b shows the spectrum of the 
noise also based on 500 samples. The overall signal-to- 
noise ratio is 0 dB. The bin signal-to-noise ratio has maxi- 
mum values of 7.5 dB in the main lobe and -6.8 dB in the 
first sidelobe. 

For a given ensemble size the spectral estimate with noise 
can be no better than the signal only spectral estimate based 
on the same number of samples. Fig. 19a shows the average 
of the signal only spectrum when M is 8, and Fig. 19b shows 
the spectral estimate when noise is added for ensemble size 
M= 8. Figures 20a and 20b are these results when M is 16. 

As expected the estimate for M = 16 is Significantly bet- 
ter than for M = 8. In addition to showing less gaps in the 
main lobe area, the M = 16 estimate has considerably less 


variation in the amplitudes of the spectral components. 
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Neither estimate shows significant components in the side- 
lobe area. This is caused by the low bin signal-to-noise 
ratio in the sidelobe area. At the peak of the sidelobe the 
Pemesignal-to-noise ratio is -6.8 dB. At this value, the 


probability of detection is 7.5 percent for M= 16. 
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V. CONCLUSIONS 


The FABLUS process provides an estimate of the spectrum 
of an unknown signal. The process performs equally well with 
narrowband and wideband signals. Two implementations of the 
FABLUS process are evaluated. One method requires twice the 
memory capacity and additional processing time, but it pro- 
duces marginally better probability of detection, and it 
requires the least restrictive assumptions about the signal 
spectrum. The accuracy of the estimate of the signal spec- 
trum is limited only by the ensemble size and the inaccuracies 


inherent in the FFT process. 
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APPENDIX A 


Tp eoeCh ee LaeOURLER TRANSFORM 


1. Equations 

A basic Fourier analysis result is that a periodic time 
function x(t) can be represented by a sum of sinusoids of 
proper frequencies and amplitudes. This sum is the Fourier 


series [Ref. 6]. 


fait) = Ay + eee + B sin(2mnt/T) ) 


N 


where T is the period of x(t). If x(t) is not periodic (or 
T -> »), then its frequency description becomes the Fourier 
transform represented mathematically as [Ref. 7]: 


F(x(t)) = xX(f) = f x(t)ese7tt 


— © 


at 


The discrete Fourier transform (DFT) uses sample values 
of x(t) obtained over a finite length of time (record length). 
The bandlimited, continuous time function x(t) having Fourier 
transform X(f) is first sampled at intervals of T seconds to 
produce the sampled time function x(mT) where m is the set 
of integers. Now, let a finite record of a sampled time func- 
tion be x(iT), i = 0,...,N-l. The Discrete Fourier transform 
of x(1iT) is 


X(q) = }x(iT) [cos(ia) -jsin (ia) ] 
at 
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where a = 2nq/N, j = v-l, q = 0O,...,N-l. Here, iisa 
subset of m, the set of integers. Dropping the T for con- 


venience and rearranging gives 


X(q) = ){x(i)cos(ia) - jx(i)sin(ia)] 
1 


which has real and imaginary parts 


Re[X(q)] = }x(i)cos(ia) 
a 
Im[X(q)] = )x(i)sin(ia). 


zx 


The squared magnitude of the transform which is used in this 


report is 


[Yx(i)cos(ia)]* + [Jx(i)sin(ia) ]? 
i i 


Ix(q) |? 





= )}x(i)x(j)cos(ia)cos(ja) + 
1) 


))*(i)x(5) sin (ia) sin(Jja) 
ij 


= J)x(i)x(3) [cos (ia) cos(ja) + 
ij 


sin(ia)sin(ja) ] 


22 ESsors 

From the sampling theorem [Ref. 8], x(mT) contains all the 
information necessary to reconstruct x(t) if T is selected 
so that 1/T is at least twice the highest non-zero frequency 


component of X(f). If X(f) exists for f > 1/2T, aliasing 
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occurs which means the frequency components of x(mT) are not 
the same as those of x(t). This is the first source of error 
faethe DFT. 

Since x(mT) contains an infinite number of values, for 
computational purposes we must truncate x(mT) to some reason- 
able length. This truncation is accomplished by selecting 
some number N of samples to be considered. Therefore, the 
record length To is NT. The truncation process is also 
called windowing because it is equivalent to multiplying x(mT) 
by a window function. In the frequency domain the multipli- 
cation of these two time functions is equivalent to convolu- 
tion of the spectrum of x(mT) with the spectrum of the window 
function. Therefore, the spectrum of x(mT) is altered by the 
window function. This is the second source of error in the 
DFT process. 

To obtain the spectrum of x(mT), the DFT assumes the avail- 
able record repeats and calculates the Fourier transform 
using the sample values. This results in values for N terms 
at frequencies which are multiples of 1/T) = 1/NT = fo: This 
discrete spectrum X (mf) may mask the maximum and minimum 
values of X(f), the true spectrum of x(t). 

The DFT output, then, consists of N complex numbers repre- 
senting the magnitudes and phases of the transform at the N 
frequencies spaced 1/NT Hz apart. However, only N/2 values 
of the transform are unique. The values from N/2 to N repre- 


sent negative frequency components [Ref. 9]. Viewed from 
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another perspective, it is unreasonable to expect a linear 
process like the DFT to produce 2N linearly independent out- 
puts (real and imaginary parts at N frequencies) with only 
N linearly independent inputs. 

We can summarize these various effects by comparing the 
DFT output with a sampled version of X(f). 

Let X" (mf) = the sampled version of X(f) 

and let X (mf )) = the DPT o@tput. 
Then, X (mf) = X" (mE ,) ~ Al(mf)) ~ Tx (mf) where Al(mf)) 
represents the distortion in X (mf) caused by aliasing and 
Tr (mf )) represents the distortion in X (mf) caused by trunca-~ 
tion. The aliasing and truncation effects can cause X (mf) 
to increase or decrease. 

The FFT is a computationally efficient algorithm used to 


generate the DFT values. 
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SEPENDEX B 
DERIVATION OF STATISTICAL PROPERTIES OF SINUSOIDAL SIGNAL 
AND UNCORRELATED, GAUSSIAN NOISE 
In this appendix, we derive the equations for the mean 
and variance of the frequency bin squared magnitudes which 
are computed by the FABLUS process. The derivation considers 
the cases of noise only and sinusoid plus noise. That is, we 
derive here E[|G.(p) |*1, var (|G. (p)|*1, E(|G.(p) [71 and 


var (|G. (p) |*1. 


1. Noise Only 

From Appendix A, if we let x(i) = n(i) = noise only and 
let Gy represent the transform under noise only conditions, 
then 


|G, (p) |* = JJn(i)n(j) [cos (ia)cos(ja) +sin(ia) sin(ja) ] 
ij 
The expected value of the squared magnitude becomes: 
E{|G.(p)|*} = JJn(i)n(j) [eos (ia) cos (ja) +sin(ia)sin(ja)] 
ahs 


If n(i) is independent from sample to sample and recognizing 
that the only random quantities are the n(i), then the ex- 
pected value is really the result of an ensemble average which 


becomes (assuming stationary noise). 


n(ijn(j) = O £OLr 1 F j 
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yn (ijn (dL) [cos* (ia) Ee Seier cies] 
i 


2 
E{ |G (p) |“3 





= Nn(i)° 


= Nn* 


Now, the squared value of the magnitude squared is 


[}X(p) 717 = { Fx (i) x (5) [eos (ia) cos (ja) tsin (ia) sin (ja) ]}° 
1) 


I 


y+} }x (i) ¥(4) x(k) x(1) [cos (ia) cos (ja) + 
i3kl 


Sin(ia)sin(ja) ] [cos(ka) cos (la) + 
sin(ka) sin(la) J 


If x(i) = n(i) = noise only, then the mean square value 


be Cenies : 
E{ (|G, (p) |*1°} = })))n(i)n (3) n(k)n(1) [cos (ia) cos (ja) + 
1 


Sin(ia)sin(ja) ] [cos (ka) cos (la) + 


sin (ka) sin(la) } 


If n(i) is independent from sample to sample and again recog- 
nizing that the n(i) are the only random quantities, then 


this expression has values only when i= j= k= 1ori=j 


and k = 1 ori=k andj2=t1orie=l1landje=k. 
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There will be N conditions when i= j = k = l. 


mific. (p)|*173 = Jn (i) *[cos* (ia) +sin* (ia) ] [cos* (ia) +sin? (ia) ] 
Z 


+ y¥n(j)n(j) nk) NK) [cos* (ja) +sin* (ja) ] [cos? (ka) +sin“ (ka) ] 
jk j#k 


+ yyn(j)N(k) (J) N(K) [cos (ja) cos (ka) +sin (ja) sin(ka) ]7 
jk jk 


+ ))ncj)n(k)n(k)n() [cos (ja) cos (ka) +sin (ja) sin(ka) ]° 
jk j#k 


The cosine squared plus sine squared terms are unity and the 


use of the produce trionometeric identity yields: 





B{{|G.(p) [717 = Jn(i)* + [JmGyatgyatkyA(K) 
i jk XK 


+ J}{tmejyntkyn (NK) ] + [n(j)n(k)n(k)n(9) 1} {cos* [ (j-k) al} 





jk i#k 
= JYn(i)* + JJn(5)2n(k)? + 
i a7 
23 Jn (4) 2n(k) “cos [ (j-k) al 
Jk jk 
E{ (|G, (p) |71°} = Nn ee Seo 2(n*) *SJcos*[ (j-k) al 


jk 


= Nn’ + (2N* -3N) (n*)? 


Now, Var(x) = eee - IE (2) 1°. Therefore, 


NHecm(ONen=anymmma) > = No(n-)- 


var (|G, (p) |71 


= nn? a (Nn? Seine) 
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If n(i) 1s Gaussian distributed then n* = ie and so 
2 2) ; 
Var(|G, (p) | } = Nn? (n?)? which is a result used to define a 


Signal-to-noise ratio in Chapter III. 


2. Signal Plus Noise 
If x(i) = s(i) + n(i) = Signal plus Noise, then letting 
Ga represent the transform under signal plus noise conditions 


the expected value of the magnitude squared becomes: 


E[|Go(p)|“] = J] Ts()*nti)1tst3) +n(3) 1 (eos (ia) cos (ja) 
ij 
+Sin(ia)sin(ij)] 
- 1 ete er oeeae revere) | Com tae) 
J 
cos (ja)+sin(ia)sin(ij)] 
If s(i) and n(i) are independent, n(i) = 0 and recognizing 


that the sine and cosine terms are not random gives 


E(|G.(p) |71 = ))}s(i)s(Jj) [cos (ia) cos (ja) +sin(ia)sin(ja) ] 
1) 


+ ))n(i)n(3j) (cos (ia) cos (ja)+sin(ia)sin(ja) ] 
13 


From previous work, the second term equals Nn? (i). ot 
s(i) = Acos(21nfci/N), then s(i) is deterministic, and we have 


by letting F = 2nfc/N 


EL|G,(p) |] = A“) Ycos (Fi) cos (F3) [cos (ia) cos (ja) +sin (ia) sin(Ja)] 
13 


+ nn 
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When F #7 p this reduces to the noise only case. 


For F = p we have 
E(|GQ(p)|*] = a*n?/4 + Nn? 


Now consider the value of the squared magnitude squared: 


Mp) -|)7 = 2b Lix(4) x (Jak) x(1) [cos (1a) cos (ja) +sin(ia) sin (ja) } 
J 
[cos (ka) cos (la) +sin(ka) sin (la) ] 
For x(i) = s(i) + n(i) we get: 
emp) |“)* = ete Ue) me Mts 221) 


[cos (ia) cos (ja)+sin(ia) sin(ja) ] 
[cos (ka) cos (la) +sin(ka)sin(la) ] 
Now, the expected value of the squared magnitude squared becomes: 


E@G.(p)|“)* = ip, =e + S(i)s(k)n(j)n(1) + 
n(i)s(j)s(k)s(1) + n(i)n(j)s(k)s(1) + 
5(i)5(j)n(k)S(1) + Sti)n(jyn(k)sS() + 
n(iys(jyn(k)Stl) + HUAN N(K)SU) + 
S(i)s(j)s(k)n(l) + S(i)s(kK)n(j)HUL) + 
n(i)S()Stk)n(1) + H(d)n(j)S(k)N(L) + 
S(i)s(jyn(kyn(l) + S(iyn(j)n(k)S(1) + 


{[cos(ia)cos(ja) + sin(ia)sin(ja) ] 


[cos (ka)cos(la) + sin(ka)sin(la) ]} 


eZ 








If n(i) has zero mean and zero third moment, and if n(i) is 
independent of s(i), and if s(i) is deterministic, and recog- 


nizing that only n(i) is random, then 


a2 


EX(/Go(p){°1°} = ))d)IsG)s(j)stk)s(Q1) + nGiynGystkys(1) + 
Take 
s(i)n(j)n(k)s(1) + n(i)s(j)n(k)s(1) + 
s(i)s(k)n(j)n(1) + n(i)s(3)s(k)n(1) + 
s(i)s(j)n(k)n(1) + n(i)n(3)n(k)n(1)] 
[cos(ia)cos(ja) + sin(ia)sin(ja) ] 
[cos(ka)cos(la) + sin(ka)sin(la) ] 
Assuming n(i)n(j) = 0 for i # j (sample to sample independence), 
then 
E{{|G.(p)|*17} = JD} fe (08 (4) (k)s (1) [cos (ia) cos (3a) + 


sin(ia)sin(ja)]{[cos(ka)cos(la) + 


sin (ka) sin(la) ] 





+ 2) Yn(4)7s (5) s (k) [eos* (ia) + sin’ (ia) ] 
ijk 


[cos (ja) cos (ka) +sin (ja) sin (ka) ] 





+ 4)J)n(i)7s(4)s(k) [cos (ia) cos (ja) + 
ijk 


sin(ia)sin(ja) [cos(ka)cos(ia) + 


sin (ka) sin(ia) ] 


+ }}})n(i)n (3) n(k)n(1) [cos(ia)cos(ja) + 
ioe 
sin(ia)sin(ja)][cos(ka)cos(la) + 


sin (ka) sin(la) ] 


13 








Assuming s(i) = Acos(Fi) where F = 2nfc/N and n(i) is Gaussian, 


then 


Bye Pe ioe EN 17 eee aiiien a 8) 


tS. (Pp) | 


~ Duane 


Pei omer) icy oer nat nie) 


Now, Var(x) = Eee ~ aeons Therefore, 


var (|G. (p) |*17 = a2nen2/2 + nN? (n7) ? 
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APPENDIX C 


PROBABILITY OF DETECTION FOR SINUSOIDAL SIGNAL AND 
UNCORRELATED, GAUSSIAN NOISE 


Using the assumptions for Method I in Chapter III, the 


probability of detection is 
Pd = P(r >Cv) 


where r is the decision variable and Cv is the critical value. 
From Appendix B for the case of white, Gaussian noise and 


Sinusoidal signal, 





= Naren 





= a°N*/4 “- nie 
var (|G,|*) = ACN? n* /2 PNG Ore 


Therefore, the decision variable become 





2 
Il? - Ic 


ye) 
it 


a°N*/4 + Nn’ - ie 


aA‘N*/4 


and 


12 





s_° = (A“N?n? /2 + N* (n°) ¢ + Né (n7} 7) /w 
SCAN ie eee nee 
5, 5 Aa noeatie tn) 17 ony. 


The critical value is: 


Cv = Z So 


z\/ (ACN 24 4n*(n*) 7) 2M) . 


Using the above expressions, 


Pd = P(r|v-u =A*N7/4) 


= Q((Cv = R)/8Q) 


= Q(z, - (A“N* YM) / (4 WV orea2n272 + 2N-{(n*)°))). 
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